Assisted crystal growing by tempering metastable vapor-liquid fluids 
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The metastable vapor-liquid coexistence of short-range attractive fluids hinders the formation 
of crystal nuclei, which in turn makes difficult the progress of the system towards its vapor-solid 
ground state. In this letter we show that crystal growth can be assisted by imposing temperature 
fluctuations. By so doing the obtained solid is nearly a fee monocrystal in contrast with the extreme 
polycrystalline structure obtained at low temperatures. The study is carried out by combining the 
replica exchange Monte Carlo method and the standard slab technique. The obtained results suggest 
a pathway for growing coherent crystals from the metastable liquid. This is particularly relevant 
for the crystallization of globular proteins. 



I. INTRODUCTION 

The temperature-density phase diagram drastically 
changes with the range of the attractive part of the in- 
terparticle potential |l|-|3(. So that, when it is less than 
half of the particle's radius, the critical point locates be- 
low the freezing curve and the corresponding vapor-liquid 
phases are metastable [il-flij]. This inhibits crystal for- 
mation. To overcome this obstacle, it was proposed to 
take advantage of natural fluctuations occurring in the 
vicinity of the critical point, which may produce the first 
stable seed to grow the crystal phase [15;] . This pathway 
for crystallization was unsuccessfully applied to globular 
proteins solutions, which are known to follow the phase 
behavior described by short-range attractive model flu- 
ids [H, 0, E3, EH]- Alternatively, temporal variations of 
protein concentration enhance crystallization, pointing 
out that fluctuations, though externally imposed, indeed 
promote nucleation and the crystal growth [16j . 

Computer simulations of short-range interacting flu- 
ids easily access and remain in the metastable region of 
the phase space, which signals a free energy barrier be- 
tween the vapor-liquid and vapor-solid states [lj], EH E3 ■ 
At sufficiently low temperatures, a liquid droplet evolves 
towards an imperfect solid made up of very small crys- 
tallites glued trough defective grain boundaries On 
the other hand, it has been shown that crystallization is 
enhanced for the patchy square well potential fluid. In 
that case crystallization takes place after the formation 
of a high density droplet of self-assembled particles and 
well outside the metastable vapor-liquid region [l9j . 

In this work we consider assisting the solid growth of 
an isotropic short-range interacting fluid by tempering 
its metastable vapor-liquid phases. In this way we ex- 
pect the system to approach to the vapor-solid ground 
state. We also study the possibility of generating the 
solid seed at the vapor phase. To deal with these issues, 



1.2 
1.0 
0.8 
10.6 
0.4 
0.2 















jl 


, , i , , , 





* godriozo@imp.mx 
t fangeles@imp.mx 
•f porea@imp.mx 



FIG. 1. Density profiles of the SW fluid for A = 1.15 and 
for different temperatures. From bottom to top in the dense 
phase, temperature geometrically decreases from T* = 0.562 
to 0.515 by a factor of 0.992. 



we employed a simulation algorithm that imposes sudden 
temperature changes over several system replicas. 



II. NUMERICAL METHOD 

Monte Carlo (MC) techniques are designed to sam- 
ple form equilibrium. However, standard MC in some 
cases cannot drive the system to its ground state due to 
the existence of deep free-energy minima. In such cases, 
the replica exchange Monte Carlo (REMC) method fre- 
quently improves the sampling [20h24| | . In this work we 
combine the REMC method with the slab technique [24[ , 
i.e., a system initially containing a liquid slab surrounded 
by its vapor is replicated for different temperatures (par- 
allel tempering) . The general idea is to simulate M repli- 
cas of the original system, being each replica at a different 
temperature, so that, the exchange of microstates among 
the cells is allowed (swap moves). These swap moves will 
make the replicas to follow some particular temperature 
routes, which can be then replicated by other simulation 
techniques. Thus, we use REMC to expose the replicas 
in the metastable state to temperature fluctuations. By 
so doing, we expect some of these replicas to escape from 
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metastability. 

Let us consider a model fluid made up of spherical 
particles at number density p, interacting through the 
attractive square- well (SW) potential, given by u(r) = oo 
for r < <7, u(r) = —e for a < r < A, and for r > A; 
being tr the particles diameter, e and A = 1.15c the well 
depth and width, respectively. For simplicity, we take 
the particles diameter as the unit length and e is given 
in units of fcgT, being ks the Boltzmann constant and 
T the absolute temperature. The corresponding reduced 
number density and temperature are given by p* = pa 3 
and T* = kgT/e, respectively. 



We employed M = 12 parallelepipeds with sides L x = 
L y = 8.0cr and L z = 40. Oct for simulating the vapor-liquid 
interface of the SW fluid. Thus, each simulation box is 
initially set with all particles randomly placed within the 
slab, in such a way that the particles concentration is 
p* ~ 0.7 and the slab is initially surrounded by vacuum. 
This initial configuration is close to a metastable liquid- 
vapor state where liquid and vapor have similar volumes, 
which allows to search for seeding at both phases. In ad- 
dition, we also set the bulk dense phase having a thick- 
ness larger than 10cr. This is, in our experience, enough 
to avoid a shift in the phase diagram due to the pres- 
ence of the interfaces between the dense and the vapor 
phases. Periodic boundary conditions are set in the three 
directions. The center of mass is placed at the box cen- 
ter. Verlet lists are implemented to improve performance. 
The highest temperature is set near but below the crit- 
ical temperature. Other temperatures are fixed by fol- 
lowing a geometrically decreasing trend. The swap ac- 
ceptance rate is, in general, close to 20 %. The replicas 
perform 10 12 trials, while maximum displacements are 
varied to yield acceptance rates close to 30 %, Long 
displacement trials are also considered. These displace- 
ments are important since they allow large jumps in the 
vapor phase with relatively large acceptance rates, while 
naturally carry out particle transference between both 
phases. After that, maximum displacements turn fixed 
and the thermodynamic properties are calculated over 
additional 10 12 configurations. 

This study consists of computing the density pro- 
files, radial distribution functions, average number of 
neighbors, and the order parameter Q§, as a func- 
tion of the temperature. The average number of first 
neighbors, N n , is the average number of pairs having 
a center-center distance smaller than 1.2<r (the vectors 
joining the centers of these pairs are named bonds). 
The order parameter Q 6 is defined as (25T - |27| Q e = 

1/2 



( III v^ m = 6 
^ 13 Z^m=- 



6m I 



where < Y 6m (9, 



> 



is the average over all bonds and configurations of the 
spherical harmonics of the angles 9 and 4> (these are the 
spherical angles of the bonds measured with respect to 
any fixed coordinate system, since Qq is invariant). 
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FIG. 2. Phase diagram obtained from the condensed and 
vapor densities of Fig. [T] (open circles). Triangles and squares 
(red and green online) are metastable vapor-liquid data taken 
from ref. [12] and ref. [13], respectively. Diamonds (purple 
online) correspond to the stable vapor-solid coexistence from 
ref. [L| (the dotted line is drown to guide the eye). The 
horizontal solid line divides the high and low temperature 
classes. 



III. RESULTS AND DISCUSSION 

From Fig. Q] it is observed that, as temperature de- 
creases, at T* « 0.533 the density profiles display an 
abrupt increase at the condensed phase while at the va- 
por phase a decrease occurs. The condensed phase of 
the fluid, with density p* d , is located at the box central 
region whereas the vapor phase, with density p v , is con- 
tiguous to both sides of the condensed phase. The profile 
corresponding to the highest temperature (T* = 0.562) 
produces the lowest value of p* d and the highest for p*, 
whereas the one corresponding to the lowest temperature 
(T* = 0.515) yields the largest p* d and smallest p*. All 
the other density profiles are for intermediate tempera- 
tures. 

Thus, the density profiles appear divided in two classes: 
The high temperature class groups the density profiles 
with p* d and p* in agreement with typical densities of liq- 
uid and vapor, respectively. This kind of density profiles 
corresponds to the metastable vapor-liquid coexistence. 
The low temperature class contains the higher values of 
p* d and lower values of p* . A distinction between the two 
classes of density profiles is the sharply defined interface 
for those at low temperature and the diffuse interface for 
those at high temperature. The latter is typical for a 
vapor-liquid coexistence, where the interface width in- 
creases with temperature. 

For the first simulation steps all density profiles behave 
like the high temperature class, since initial configura- 
tions are set in metastable states, i. e., a liquid slab with 
density close to 0.7 and its corresponding vapor. At this 
stage, all swap trials have acceptance rates over 20%. 
After a relatively large number of steps, driven by the 
temperature variations, a given replica escapes from its 
initial state towards its ground state, which give rise to 
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the second kind of profiles. The breaking away from the 
metastable state is accompanied by a strong decrease of 
the configuration energy, which makes the REMC tech- 
nique to preferable locate it at the lowest temperature. 
A similar process occurs for other replicas until the re- 
maining high temperature replicas unlikely access tem- 
peratures below T* = T* = 0.533 ± 0.002. In parallel, 
a free energy gap among the replicas at the metastable 
region and those closer to their ground state is devel- 
oped. Due to this gap the swap moves between the high 
and low temperature replicas turn practically forbidden. 
Consequently, the high temperature replicas do not ac- 
cess low temperatures states and so, they cannot escape 
from the vapor-liquid metastable region. Under this ra- 
tionale, not all replicas break away towards the ground 
state, but just some of them, implying that the whole 
system of replicas does not reach equilibrium. However, 
after approximately 3 x 10 11 trials the whole system be- 
comes stationary and thus, sampling is carried out. 

The well defined condensed-phase and vapor regions 
directly lead to the vapor and condensed densities, p* 
and p* d , by taking average at the different regions of the 
profiles. This is done without the need of fitting a hy- 
perbolic function, as it is commonly used. Fig. [2] shows 
the coexistence phase diagram obtained by plotting p* 
and p* d as a function of temperature. The horizontal line 
points out the discontinuity of the p* and p* d branches 
at T* w 0.533. Such a discontinuity separates the two 
density profile classes discussed in Fig. [TJ 

For the upper part of the diagram we employed the 
exponential fitting (with j3 — 0.325) to estimate the crit- 
ical point [23. We obtained T* = 0.575 and p* = 0.435 
for the critical temperature and density, respectively. As 
temperature decreases (in Fig. [5]), there is an abrupt shift 
from vapor-liquid to a denser phase-vapor coexistence. 
For T* < 0.533 we get the p* d and p* results shown under 
the horizontal line. These branches correspond to the 
low temperature class concentration profiles in Fig. [1] 
The density values of these branches indicate that the 
low temperature condensed phase is a solid. This means 
that the imposed temperature drops on the replicas aid 
the system to avoid the free energy barrier between the 
metastable and the ground state. That is, the crystal 
seeds are produced at low temperatures, which makes 
some of these microstates to appear at the other side of 
the free energy barrier when they are reset at high tem- 
perature in the swapping process. 

Above the discontinuity, our results well agree with 
those previously reported and carried out without inter- 
faces [HI, [H[ . The results of ref . 12 1 are shown as tri- 
angles in Fig. [2] Our results also agree with previous 
calculations using the slab technique [17|. Below the dis- 
continuity, our data approach to the stable vapor-solid 
region described by Liu et al. (diamonds in Fig. [5]). Ac- 
tually, the vapor density branch falls over the interpola- 
tion of Liu's data, which were obtained by imposing a 
vapor-solid interface [HI, [29| . In our case the vapor-solid 
interface is spontaneously formed. The observed agree- 
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FIG. 3. The first neighbor number, N n , in the condensed 
region as a function of temperature (top). Order parameter 
Qg as a function of temperature (bottom). The vertical solid 
line corresponds to T a *. 



ment and the free energy decreasing confirms our asser- 
tion that above and below T* the data correspond to the 
metastable and stable coexistence diagrams, respectively. 

To characterize the high density phase, the first neigh- 
bor number, N n , and order parameter Qq were followed 
as a function of temperature. The results are shown in 
Fig. [31 where a vertical line is plotted at T*. These 
quantities were obtained only for the condensed phase, 
and avoiding the interfaces. For T* > T*, we obtained 
N n « 7, and Qq < 0.1, which are the typical values 
for a liquid. At T*, N n and Qq display a discontinuity: 
N n jumps from 7.95 to 11.20, whereas Qq from 0.064 to 
0.400. It is kno wn that the order parameter, Q§, tends 
to 1/ y f N s N n /2 for a completely random arrangement of 
particles, 0.5745 for a face cubic centered (fee) crystal, 
0.4848 for a hexagonal close packed (hep) structure, and 
0.3536 for a simple cubic geometry (being N s the number 
of pairs considered for the calculation) . For T* < T* we 
get 0.40 <Q 6 < 0.51 and 11.2 < N n < 12.0, which signa- 
tures the presence of a solid condensed phase constituted 
by fee or/and hep crystal domains. 

We gain further insight of the condensed phase struc- 
ture by building the radial distribution functions (RDFs). 
This is shown in Fig. [4] for the lowest and highest com- 
puted temperatures. The high temperature RDF shows 
the typical behavior for a SW liquid, i. e., there are two 
discontinuities at r — a and r = a + A, followed by two 
broad peaks at rja k, 1.85 and rja « 2.06. At the lo- 
cation of the SW the RDF has large values, being the 
largest at contact, r — a, and the lowest at r = a + A. 
There is a nearly linear decreasing trend at the SW re- 
gion. On the other hand, the RDF for the lowest tem- 
perature shows several well defined peaks, pointing out 
a long range order. In this case, the highest peak is at 
rja = 1.06, and the following are at rja = 1.51, 1.85, 
2.13, and 2.40, respectively, which are slightly shifted to 
the right with respect to those corresponding to a fee 
crystal [2^. This is in common for all obtained solid 
phases. This result contrasts with the rapid decay of 



4 







•51 






o 


^ k 




JLwA } 





r 



FIG. 4. Radial distribution functions (RDFs) for the con- 
densed phases. The dark and light (blue) lines correspond to 
the lowest and highest computed temperatures, respectively. 
The insets are the corresponding snapshots for the vapor- 
liquid (top) and the vapor-solid (bottom) interfaces. 



the vapor-liquid (top), and vapor-solid (bottom) coex- 
isting phases in the inset of Fig. |4] In these snapshots 
many of the effects discussed throughout the paper are 
visualized. That is, the sharply defined vapor-solid inter- 
face (bottom) , the fee solid structure (bottom) , the liquid 
structure (top), and a more dense vapor phase coexist- 
ing with the liquid than with the solid (top and bottom) . 
The crystalline structure explains the presence and ab- 
sence of oscillations in the vapor-solid density profiles of 
Fig. [I] If any of the crystal planes is parallel to the inter- 
face (xy-plane), oscillations appear, whereas they vanish 
in any other case. This is why several low temperature 
profiles in Fig. [1] are flat even though they correspond 
to crystalline structures. The crystal growth cannot be 
controlled by the method, and so, the crystal planes may 
or may not be parallel to the interface. The snapshots 
also show the aggregation of vapor particles at both tem- 
peratures. However, no crystal nuclei are observed in the 
vapor phase in any case though the vapor is set relatively 
large in all cases. 



the RDF for the extreme polycrystalline solid produced 
in the absence of temperature fluctuations, which has a 
density slightly above that of the liquid and a fractal 
structure [18| . 

In practice, an equivalent procedure of the REMC sim- 
ulation technique could be reached by successively im- 
posing increasing and decreasing variations of temper- 
ature on a colloidal suspension. So that, in this way 
one could expect to produce a crystal seed at low tem- 
perature, which then can grow with the imposed fluc- 
tuating temperature. The growing process incorporates 
particles and new grains at the boundary of the grow- 
ing nucleus [l8j]- At this point, a large temperature aids 
reorienting the adjacent grains and relaxing their bound- 
aries, however, a large exposure may lead to its total 
melting. On the other hand, a low temperature fastens 
the growing processes with the consequent incorporation 
of defects. Hence, tempering combines the reorienting, 
boundary relaxing, and growing processes of the nuclei, 
yielding a more coherent solid. Temperature fluctuations 
may also aid with subsequent grain coalescence processes. 

Finally, we included two representative snapshots of 



IV. CONCLUSIONS 

We studied the response of the SW fluid metastable 
liquid-vapor coexistence to temperature fluctuations un- 
der the scheme of the replica exchange Monte Carlo al- 
gorithm. Even though we imposed a large vapor volume, 
the solid seeding was observed always inside the liquid 
phase. This supports the hypothesis that the formation 
of a dense droplet is a necessary condition to produce the 
crystal 19]. The obtained solid is in all cases a nearly fee 
monocrystal in opposition to the highly polycrystalline 
structure obtained by a simple quenching process fl8j . 
Our results suggest that following the seed formation at 
low temperatures, the crystal phase can grow and relax 
approaching to the ground state when temperature is in- 
creased and decreased successively. 
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